Eggshell_review_project
Lucas Legendre
14/12/2021
Compiled under R version 4.1.2 (2021-11-01)
WARNING: edit the working directory to your preferred folder.
This document details all analyses performed in R for the study:
Legendre, L. J., S. Choi, and J. A. Clarke. The diverse terminology and complex evolution of reptile eggshell microstructure.
For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.
Loading packages
library(ape)
library(castor)
library(evobiR)
library(ggtree)
library(phytools)
library(RColorBrewer)Data and tree
- Tree
tree<-read.nexus("treewhole_newversion.trees.nex")
plotTree(tree, fsize=0.5,lwd=1,type="fan",ftype="i")
- Data
data<-read.table("Datawhole_newproject.txt", header=T)
setdiff(tree$tip.label, data$Taxon) # taxa and data match## character(0)rownames(data)<-data$Taxon
datanew<-ReorderData(tree, data)Ancestral state reconstruction: discrete trait (soft/semi-rigid/hard)
- Extract data vectors for each coding
thisstudy<-datanew$Type_2021; names(thisstudy)<-rownames(datanew)
norell<-datanew$Type_Norell; names(norell)<-rownames(datanew)
legendre<-datanew$Type_Legendre; names(legendre)<-rownames(datanew)
# Color palette
cols<-setNames(c("royalblue","green3","red3"),c("Soft","Semi-rigid","Hard"))
# Visualize tree with node numbers
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
Using character coding defined in this study
- Perform SIMMAP using
phytools– 1000 iterations, using AIC to select best model (code modified from Liam Revell, see here)
x<-thisstudy
aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
d.aic<-aic-min(aic)
exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010672548 0.0005336274 0.0005336274
## Semi-rigid 0.0005336274 -0.0010672548 0.0005336274
## Soft 0.0005336274 0.0005336274 -0.0010672548
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010457689 0.0003790695 0.0006666995
## Semi-rigid 0.0003790695 -0.0008898238 0.0005107544
## Soft 0.0006666995 0.0005107544 -0.0011774539
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.00000000 0.0000000000
## Semi-rigid 0.0058104610 -0.00962594 0.0038154791
## Soft 0.0004728781 0.00000000 -0.0004728781
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.logL## ER SYM ARD
## -56.06148 -55.75401 -43.99790AIC<-mapply(aic,logL,c(1,3,6))
AIC## ER SYM ARD
## 114.1230 117.5080 99.9958AIC.W<-aic.w(AIC)
AIC.W## ER SYM ARD
## 0.0008548391 0.0001573372 0.9989878237nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 1 0 999o1<-make.simmap(tree,x,model="ARD",nsim=999)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.00000000 0.0000000000
## Semi-rigid 0.0058104610 -0.00962594 0.0038154791
## Soft 0.0004728781 0.00000000 -0.0004728781
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.o2<-make.simmap(tree,x,model="ER",nsim=1)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0010672548 0.0005336274 0.0005336274
## Semi-rigid 0.0005336274 -0.0010672548 0.0005336274
## Soft 0.0005336274 0.0005336274 -0.0010672548
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treesstudy<-c(o1,o2)
objstudy<-describe.simmap(treesstudy)
objstudy## 1000 trees with a mapped discrete character with states:
## Hard, Semi-rigid, Soft
##
## trees have 18.211 changes between states on average
##
## changes are of the following types:
## Hard,Semi-rigid Hard,Soft Semi-rigid,Hard Semi-rigid,Soft Soft,Hard
## x->y 0.003 0.001 10.112 6.631 1.462
## Soft,Semi-rigid
## x->y 0.002
##
## mean total time spent in each state is:
## Hard Semi-rigid Soft total
## raw 8019.4557598 1739.0045051 3095.1775697 12853.64
## prop 0.6239055 0.1352928 0.2408017 1.00plot(objstudy,type="fan",fsize=0.01,lwd=1,ftype="i", colors=cols,ylim=c(-2,Ntip(tree)),offset=20, part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Here, a semi-rigid eggshell is the ancestral state for almost all major reptile clades – reptiles, lepidosaurs, squamates, turtles, archelosaurs, archosaurs, dinosaurs,ornithischians, and saurischians. However, there are only 6 taxa with semi-rigid eggshells out of 203!
=> Some of these taxa may have a excessive influence on this result – probably the two sauropodomorphs with that eggshell type, since they are the closer in age to the root and to the aforementioned subclades.
Using character coding from Norell et al. (2020)
x<-norell
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013811074 0.0006905537 0.0006905537
## Semi-rigid 0.0006905537 -0.0013811074 0.0006905537
## Soft 0.0006905537 0.0006905537 -0.0013811074
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013651414 0.0004938598 0.0008712815
## Semi-rigid 0.0004938598 -0.0010970016 0.0006031418
## Soft 0.0008712815 0.0006031418 -0.0014744233
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.000000000 0.000000000 0.000000000
## Semi-rigid 0.007423526 -0.012130029 0.004706502
## Soft 0.001380656 0.001222171 -0.002602827
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.logL## ER SYM ARD
## -68.56111 -68.13921 -55.04396AIC<-mapply(aic,logL,c(1,3,6))
AIC## ER SYM ARD
## 139.1222 142.2784 122.0879AIC.W<-aic.w(AIC)
AIC.W## ER SYM ARD
## 1.999592e-04 4.126532e-05 9.997588e-01nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 0 0 1000treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.000000000 0.000000000 0.000000000
## Semi-rigid 0.007423526 -0.012130029 0.004706502
## Soft 0.001380656 0.001222171 -0.002602827
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.objnorell<-describe.simmap(treesnorell)
objnorell## 1000 trees with a mapped discrete character with states:
## Hard, Semi-rigid, Soft
##
## trees have 23.847 changes between states on average
##
## changes are of the following types:
## Hard,Semi-rigid Hard,Soft Semi-rigid,Hard Semi-rigid,Soft Soft,Hard
## x->y 0 0 8.254 5.153 5.562
## Soft,Semi-rigid
## x->y 4.878
##
## mean total time spent in each state is:
## Hard Semi-rigid Soft total
## raw 7730.2217576 1.110435e+03 4012.9805883 12853.64
## prop 0.6014034 8.639076e-02 0.3122058 1.00plot(objnorell,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Completely different result – let us check which taxa are different from the previous ASR…
datanew[c(which(datanew$Type_2021!=datanew$Type_Norell)),]## Clade Taxon
## Massospondylus Non-avian_dinosaurs Massospondylus
## Lufengosaurus Non-avian_dinosaurs Lufengosaurus
## Chelonoidis_niger Chelonians Chelonoidis_niger
## Testudo_marginata Chelonians Testudo_marginata
## Testudo_graeca Chelonians Testudo_graeca
## Stigmochelys_pardalis Chelonians Stigmochelys_pardalis
## Caretta_caretta Chelonians Caretta_caretta
## Eggshell_thickness Egg_mass Type_2021 Type_Norell
## Massospondylus 100 158.00000 Semi-rigid Soft
## Lufengosaurus 85 500.00000 Semi-rigid Soft
## Chelonoidis_niger 400 101.64000 Hard Soft
## Testudo_marginata 90 16.63200 Hard Soft
## Testudo_graeca 165 16.51650 Hard Soft
## Stigmochelys_pardalis 210 6.37875 Hard Semi-rigid
## Caretta_caretta 80 37.20087 Soft Semi-rigid
## Type_Legendre
## Massospondylus Hard
## Lufengosaurus Hard
## Chelonoidis_niger Hard
## Testudo_marginata Hard
## Testudo_graeca Hard
## Stigmochelys_pardalis Hard
## Caretta_caretta Soft# Only 7 taxa are different – two non-avian dinosaurs and five turtles
# Changing character state for the two non-avian dinosaurs
x[c(21,22)]<-"Semi-rigid"
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.00144622 0.00072311 0.00072311
## Semi-rigid 0.00072311 -0.00144622 0.00072311
## Soft 0.00072311 0.00072311 -0.00144622
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0013771501 0.0005401113 0.0008370389
## Semi-rigid 0.0005401113 -0.0013929825 0.0008528712
## Soft 0.0008370389 0.0008528712 -0.0016899101
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.0000000000 0.000000000
## Semi-rigid 0.0068533111 -0.0122551684 0.005401857
## Soft 0.0006039133 0.0006682603 -0.001272174
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 0 0 1000treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard 0.0000000000 0.0000000000 0.000000000
## Semi-rigid 0.0068533111 -0.0122551684 0.005401857
## Soft 0.0006039133 0.0006682603 -0.001272174
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.objnorell<-describe.simmap(treesnorell)
plot(objnorell,fsize=0.5,lwd=1,ftype="i",colors=cols,ylim=c(-2,Ntip(tree)),part=0.95)
add.simmap.legend(colors=cols,prompt=FALSE,x=20,y=180)
All major clades are now recovered as semi-rigid again => strong influence of the two sauropodomorphs.
Using character coding from Legendre et al. (2020)
x<-legendre
logL<-sapply(c("ER","SYM","ARD"),
function(model,tree,x) make.simmap(tree,x,model)$logL,
tree=tree,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007636699 0.0003818350 0.0003818350
## Semi-rigid 0.0003818350 -0.0007636699 0.0003818350
## Soft 0.0003818350 0.0003818350 -0.0007636699
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006269773 0.0000000000 0.0006269773
## Semi-rigid 0.0000000000 -0.0005633914 0.0005633914
## Soft 0.0006269773 0.0005633914 -0.0011903687
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001519306 0.0000000000 0.0001519306
## Semi-rigid 0.0019150212 -0.0019150212 0.0000000000
## Soft 0.0019555075 0.0005401052 -0.0024956127
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.logL## ER SYM ARD
## -45.50400 -43.12910 -41.26503AIC<-mapply(aic,logL,c(1,3,6))
AIC## ER SYM ARD
## 93.00800 92.25820 94.53006AIC.W<-aic.w(AIC)
AIC.W## ER SYM ARD
## 0.3422281 0.4978884 0.1598835nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 342 498 160l1<-make.simmap(tree,x,model="ER",nsim=342)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007636699 0.0003818350 0.0003818350
## Semi-rigid 0.0003818350 -0.0007636699 0.0003818350
## Soft 0.0003818350 0.0003818350 -0.0007636699
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.l2<-make.simmap(tree,x,model="SYM",nsim=498)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006269773 0.0000000000 0.0006269773
## Semi-rigid 0.0000000000 -0.0005633914 0.0005633914
## Soft 0.0006269773 0.0005633914 -0.0011903687
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.l3<-make.simmap(tree,x,model="ARD",nsim=160)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001519306 0.0000000000 0.0001519306
## Semi-rigid 0.0019150212 -0.0019150212 0.0000000000
## Soft 0.0019555075 0.0005401052 -0.0024956127
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treeslegendre<-c(l1,l2,l3)
objlegendre<-describe.simmap(treeslegendre)
objlegendre## 1000 trees with a mapped discrete character with states:
## Hard, Semi-rigid, Soft
##
## trees have 10.214 changes between states on average
##
## changes are of the following types:
## Hard,Semi-rigid Hard,Soft Semi-rigid,Hard Semi-rigid,Soft Soft,Hard
## x->y 0.411 4.143 0.14 0.061 3.772
## Soft,Semi-rigid
## x->y 1.687
##
## mean total time spent in each state is:
## Hard Semi-rigid Soft total
## raw 9104.1361000 220.33623291 3529.1655017 12853.64
## prop 0.7082926 0.01714194 0.2745655 1.00plot(objlegendre,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Reptiles, lepidosaurs, and squamates are recovered as ancestrally soft-shelled, while turtles, archelosaurs, archosaurs, and all dinosaur clades are recovered as hard-shelled.
The pattern seems to follow whichever character state is coded for the two sauropodomorphs Lufengosaurus and Massospondylus.
To test this hypothesis, we must remove all other taxa (n = 7) that are not coded the same in all three studies, and see if the pattern is still present.
- Remove all taxa with differences in coding, except the two sauropodomorphs
remove<-rownames(datanew[c(which(datanew$Type_Legendre!=datanew$Type_Norell)),][c(3:9),])
x<-thisstudy[!names(thisstudy)%in%remove]
treenew<-drop.tip(tree, setdiff(tree$tip.label, names(x)))
# Visualize tree with node numbers
ggtree(treenew) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()
Replicate analysis on the new tree 3 times, with identical coding for all taxa except the two sauropodomorphs
- Coded as semi-rigid
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0009105572 0.0004552786 0.0004552786
## Semi-rigid 0.0004552786 -0.0009105572 0.0004552786
## Soft 0.0004552786 0.0004552786 -0.0009105572
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007204487 0.0000000000 0.0007204487
## Semi-rigid 0.0000000000 -0.0007881215 0.0007881215
## Soft 0.0007204487 0.0007881215 -0.0015085703
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001356366 0.000000000 0.0001356366
## Semi-rigid 0.0033499438 -0.003349944 0.0000000000
## Soft 0.0016634621 0.000864281 -0.0025277432
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 75 97 828s1<-make.simmap(treenew,x,model="ER",nsim=75)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0009105572 0.0004552786 0.0004552786
## Semi-rigid 0.0004552786 -0.0009105572 0.0004552786
## Soft 0.0004552786 0.0004552786 -0.0009105572
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.s2<-make.simmap(treenew,x,model="SYM",nsim=97)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007204487 0.0000000000 0.0007204487
## Semi-rigid 0.0000000000 -0.0007881215 0.0007881215
## Soft 0.0007204487 0.0007881215 -0.0015085703
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.s3<-make.simmap(treenew,x,model="ARD",nsim=828)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001356366 0.000000000 0.0001356366
## Semi-rigid 0.0033499438 -0.003349944 0.0000000000
## Soft 0.0016634621 0.000864281 -0.0025277432
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treenewsr<-c(s1,s2,s3)
objsr<-describe.simmap(treenewsr)
plot(objsr,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Most inner clades are now recovered as soft-shelled – most likely an artefact due to the position of Mussaurus – closer to the other inner nodes than any other taxon, since it is the oldest specimen in the sample.
- Coded as soft
x[c(21,22)]<-"Soft"
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008576981 0.0004288491 0.0004288491
## Semi-rigid 0.0004288491 -0.0008576981 0.0004288491
## Soft 0.0004288491 0.0004288491 -0.0008576981
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007184639 0.0000000000 0.0007184639
## Semi-rigid 0.0000000000 -0.0005102417 0.0005102417
## Soft 0.0007184639 0.0005102417 -0.0012287055
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001351033 0.0000000000 0.0001351033
## Semi-rigid 0.0024265091 -0.0024265091 0.0000000000
## Soft 0.0017286691 0.0005520113 -0.0022806804
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 47 105 848so1<-make.simmap(treenew,x,model="ER",nsim=46)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0008576981 0.0004288491 0.0004288491
## Semi-rigid 0.0004288491 -0.0008576981 0.0004288491
## Soft 0.0004288491 0.0004288491 -0.0008576981
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.so2<-make.simmap(treenew,x,model="SYM",nsim=106)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007184639 0.0000000000 0.0007184639
## Semi-rigid 0.0000000000 -0.0005102417 0.0005102417
## Soft 0.0007184639 0.0005102417 -0.0012287055
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.so3<-make.simmap(treenew,x,model="ARD",nsim=848)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001351033 0.0000000000 0.0001351033
## Semi-rigid 0.0024265091 -0.0024265091 0.0000000000
## Soft 0.0017286691 0.0005520113 -0.0022806804
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treenews<-c(so1, so2, so3)
objs<-describe.simmap(treenews)
plot(objs,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Same result as with semi-rigid coding, unsurprisingly.
- Coded as hard
x[c(21,22)]<-"Hard"
logL<-sapply(c("ER","SYM","ARD"),
function(model,treenew,x) make.simmap(treenew,x,model)$logL,
tree=treenew,x=x)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007868496 0.0003934248 0.0003934248
## Semi-rigid 0.0003934248 -0.0007868496 0.0003934248
## Soft 0.0003934248 0.0003934248 -0.0007868496
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006471696 0.000000000 0.0006471696
## Semi-rigid 0.0000000000 -0.000567363 0.0005673630
## Soft 0.0006471696 0.000567363 -0.0012145326
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001501237 0.0000000000 0.0001501237
## Semi-rigid 0.0018901522 -0.0018901522 0.0000000000
## Soft 0.0020179175 0.0005399515 -0.0025578690
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim## ER SYM ARD
## 354 496 150h1<-make.simmap(treenew,x,model="ER",nsim=354)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0007868496 0.0003934248 0.0003934248
## Semi-rigid 0.0003934248 -0.0007868496 0.0003934248
## Soft 0.0003934248 0.0003934248 -0.0007868496
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.h2<-make.simmap(treenew,x,model="SYM",nsim=496)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0006471696 0.000000000 0.0006471696
## Semi-rigid 0.0000000000 -0.000567363 0.0005673630
## Soft 0.0006471696 0.000567363 -0.0012145326
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.h3<-make.simmap(treenew,x,model="ARD",nsim=150)## make.simmap is sampling character histories conditioned on
## the transition matrix
##
## Q =
## Hard Semi-rigid Soft
## Hard -0.0001501237 0.0000000000 0.0001501237
## Semi-rigid 0.0018901522 -0.0018901522 0.0000000000
## Soft 0.0020179175 0.0005399515 -0.0025578690
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Hard Semi-rigid Soft
## 0.3333333 0.3333333 0.3333333## Done.treenewh<-c(h1, h2, h3)
objh<-describe.simmap(treenewh)
plot(objh,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)
Ancestral state for dinosaurs and archosaurs, and all major less inclusive clades except pterosaurs, is hard-shelled, as in the Legendre et al. (2020) coding.
To check how strongly are these results influenced by branch length information, we can replicate these analyses using maximum parsimony, which does not consider branch length information, using castor.
- Coding from this study
studyN<-thisstudy
studyN[studyN=="Soft"]<-1; studyN[studyN=="Semi-rigid"]<-2; studyN[studyN=="Hard"]<-3
studyN<-as.numeric(studyN); names(studyN)<-names(thisstudy)
MPS<-asr_max_parsimony(tree, studyN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPS$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(studyN,sort(unique(studyN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)
MPS2<-MPS$ancestral_likelihoods; rownames(MPS2)<-c(204:369)- Coding from Norell et al. (2020)
norellN<-norell
norellN[norellN=="Soft"]<-1; norellN[norellN=="Semi-rigid"]<-2; norellN[norellN=="Hard"]<-3
norellN<-as.numeric(norellN); names(norellN)<-names(norell)
MPN<-asr_max_parsimony(tree, norellN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPN$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(norellN,sort(unique(norellN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=180,fsize=0.8)
MPN2<-MPN$ancestral_likelihoods; rownames(MPN2)<-c(204:369)- Coding from Legendre et al. (2020)
legendreN<-legendre
legendreN[legendreN=="Soft"]<-1; legendreN[legendreN=="Semi-rigid"]<-2; legendreN[legendreN=="Hard"]<-3
legendreN<-as.numeric(legendreN); names(legendreN)<-names(legendre)
MPL<-asr_max_parsimony(tree, legendreN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPL$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(legendreN,sort(unique(legendreN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)
MPL2<-MPL$ancestral_likelihoods; rownames(MPL2)<-c(204:369)- With the reduced tree, changing only the two sauropodomorphs
xN<-x
xN[xN=="Soft"]<-1; xN[xN=="Semi-rigid"]<-2; xN[xN=="Hard"]<-3
xN<-as.numeric(xN); names(xN)<-names(x)
# Coded as soft
xN[c(21,22)]<-1
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods; rownames(MPx2)<-c(197:355)
# Coded as semi-rigid
xN[c(21,22)]<-2
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods; rownames(MPx2)<-c(197:355)
# Coded as hard
xN[c(21,22)]<-3
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)
MPx2<-MPx$ancestral_likelihoods; rownames(MPx2)<-c(197:355)Ancestral state reconstruction – continuous character (eggshell thickness)
- Data and tree
treedi<-multi2di(tree, random=FALSE)
AET<-log1p(datanew$Eggshell_thickness); names(AET)<-rownames(datanew)
RET<-log1p(datanew$Eggshell_thickness/datanew$Egg_mass); names(RET)<-rownames(datanew)- Phylogenetic signal (
contMapassumes a Brownian motion model)
phylosig(treedi, AET, method="lambda", test=TRUE)##
## Phylogenetic signal lambda : 0.999957
## logL(lambda) : -309.493
## LR(lambda=0) : 245.84
## P-value (based on LR test) : 2.09625e-55phylosig(treedi, RET, method="lambda", test=TRUE)##
## Phylogenetic signal lambda : 0.999934
## logL(lambda) : -228.063
## LR(lambda=0) : 197.053
## P-value (based on LR test) : 9.18051e-45Very strong signal for both absolute and relative eggshell thickness.
- ASR for absolute eggshell thickness
phyloplotAET<-contMap(treedi, AET, plot=F)
plot(setMap(phyloplotAET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
prompt=FALSE,x=-300,y=280)
# Ancestral states
expm1(fastAnc(treedi, AET))## Ancestral character estimates using fastAnc:
## 204 205 206 207 208 209
## 28.079153 28.079153 26.118366 19.854323 36.974371 57.309311
## 210 211 212 213 214 215
## 59.388937 38.950551 37.943244 62.261450 63.628751 60.801065
## 216 217 218 219 220 221
## 18.564828 17.370830 16.914159 21.363127 14.483029 17.067109
## 222 223 224 225 226 227
## 18.167139 27.778883 11.637074 6.394273 5.261709 11.868343
## 228 229 230 231 232 233
## 10.874468 12.260022 12.122288 14.011868 16.557668 31.172419
## 234 235 236 237 238 239
## 13.895415 14.226503 16.834523 16.374609 14.489393 8.215824
## 240 241 242 243 244 245
## 8.889975 16.110137 17.584868 17.641372 11.369380 11.391237
## 246 247 248 249 250 251
## 11.091746 7.684418 36.692394 129.109053 132.829774 140.343937
## 252 253 254 255 256 257
## 145.689626 141.878765 145.186350 61.134926 159.995521 188.575290
## 258 259 260 261 262 263
## 202.108248 218.163318 174.871933 149.316431 250.541379 333.057440
## 264 265 266 267 268 269
## 379.854236 544.032218 573.191181 610.483803 607.697486 498.358959
## 270 271 272 273 274 275
## 33.924264 333.743365 451.330862 630.020508 397.177198 462.000397
## 276 277 278 279 280 281
## 470.274154 478.924007 451.589838 28.662359 2.833103 3.998462
## 282 283 284 285 286 287
## 11.270283 3.351885 0.014325 32.566033 69.398144 453.836756
## 288 289 290 291 292 293
## 1437.961087 1437.961087 453.836756 1834.612767 1834.612767 453.836756
## 294 295 296 297 298 299
## 1523.667987 1544.970306 27.280341 13.898784 72.135564 1.763827
## 300 301 302 303 304 305
## 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241
## 306 307 308 309 310 311
## 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241 1601.759241
## 312 313 314 315 316 317
## 1601.759241 1601.759241 311.981337 473.591727 855.767272 855.767272
## 318 319 320 321 322 323
## 855.767272 531.023321 533.631804 533.631804 555.967722 1285.161706
## 324 325 326 327 328 329
## 566.272234 605.254074 605.254074 605.254074 605.254074 1101.569284
## 330 331 332 333 334 335
## 1391.646422 539.617230 539.617230 539.617230 655.621109 418.158680
## 336 337 338 339 340 341
## 418.158680 418.158680 399.605014 399.605014 399.605014 392.933851
## 342 343 344 345 346 347
## 399.605014 383.372932 482.275717 482.275717 482.275717 482.275717
## 348 349 350 351 352 353
## 482.275717 487.055112 482.275717 672.352457 688.294201 743.214121
## 354 355 356 357 358 359
## 749.815848 642.557934 1304.082510 319.361966 778.723541 908.510163
## 360 361 362 363 364 365
## 820.749499 423.101767 255.012307 300.635949 343.596583 298.819020
## 366 367 368 369 370 371
## 290.372065 222.110312 217.727290 177.850464 166.734252 161.025519
## 372 373 374 375 376 377
## 229.672210 234.226968 236.209118 233.541444 232.183846 231.144902
## 378 379 380 381 382 383
## 196.457289 216.843344 210.907445 208.757659 210.716135 206.846733
## 384 385 386 387 388 389
## 201.747492 251.211818 289.690559 316.985021 253.493514 257.535012
## 390 391 392 393 394 395
## 243.362806 151.829900 155.868914 164.019364 144.779678 128.402727
## 396 397 398 399 400 401
## 106.958920 100.296111 118.838462 109.666788 73.949794 103.666941
## 402 403 404 405
## 99.285261 98.490117 103.797286 101.497666Ancestral eggshell thickness seems to be intermediate for all major nodes, including archosaurs and dinosaurs. However, many extant lepidosaurs are also recovered as intermediate – likely due to the extremely low values in pterosaurs, reported as lacking a calcareous layer entirely, which shift any thicker eggshell towards the middle of the spectrum.
Furthermore, the pattern strongly follows that of body mass, as expected (e.g. Stein et al., 2019; Legendre and Clarke, 2021).
- ASR for relative eggshell thickness
phyloplotRET<-contMap(treedi, RET, plot=F)
plot(setMap(phyloplotRET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
prompt=FALSE,x=-300,y=280)
# Ancestral states
expm1(fastAnc(treedi, RET))## Ancestral character estimates using fastAnc:
## 204 205 206 207 208 209 210 211
## 3.409640 3.409640 4.025467 7.755305 20.354753 37.931589 39.819215 25.689681
## 212 213 214 215 216 217 218 219
## 39.657045 42.671676 44.850528 53.000367 7.500649 7.357594 7.168380 6.552197
## 220 221 222 223 224 225 226 227
## 9.462534 7.377703 1.646676 1.002991 1.610506 0.936549 0.860297 1.828548
## 228 229 230 231 232 233 234 235
## 1.596624 2.012481 3.474753 1.682958 2.085547 3.133932 1.180141 1.164824
## 236 237 238 239 240 241 242 243
## 7.649844 9.184624 11.645558 12.734058 13.645813 13.020440 16.528768 20.143413
## 244 245 246 247 248 249 250 251
## 21.874101 24.369790 23.820456 14.769006 2.599039 6.626096 6.904628 8.138380
## 252 253 254 255 256 257 258 259
## 10.904586 10.965487 12.341525 2.047360 11.403766 11.235749 11.089139 11.128995
## 260 261 262 263 264 265 266 267
## 10.345142 8.942119 11.414881 8.881152 9.028346 10.830652 10.810594 10.535716
## 268 269 270 271 272 273 274 275
## 10.641304 9.256647 2.124890 4.558317 5.340322 6.761045 4.477445 5.802632
## 276 277 278 279 280 281 282 283
## 5.664138 5.483004 6.059812 1.692954 0.457460 0.550221 1.153386 0.408094
## 284 285 286 287 288 289 290 291
## 0.003995 1.252572 1.164206 1.585486 5.318273 5.318273 1.585486 0.666241
## 292 293 294 295 296 297 298 299
## 0.666241 1.585486 1.161423 1.234126 0.894187 0.428783 0.388280 0.130108
## 300 301 302 303 304 305 306 307
## 0.766339 0.766339 0.766339 0.766339 0.766339 0.766339 0.766339 0.766339
## 308 309 310 311 312 313 314 315
## 0.766339 0.766339 0.766339 0.766339 0.766339 0.766339 2.783888 4.262094
## 316 317 318 319 320 321 322 323
## 4.206638 4.206638 4.206638 5.092815 5.314153 5.314153 5.246828 2.804377
## 324 325 326 327 328 329 330 331
## 5.263532 4.918171 4.918171 4.918171 4.918171 3.280773 2.041443 5.646772
## 332 333 334 335 336 337 338 339
## 5.646772 5.646772 4.638261 8.395933 8.395933 8.395933 8.720703 8.720703
## 340 341 342 343 344 345 346 347
## 8.720703 8.804955 8.720703 8.943290 8.062171 8.062171 8.062171 8.062171
## 348 349 350 351 352 353 354 355
## 8.062171 8.106581 8.062171 3.124877 3.017738 2.342408 2.259157 2.383106
## 356 357 358 359 360 361 362 363
## 0.661660 4.634275 2.112184 1.859717 1.878840 1.176578 10.833970 9.311577
## 364 365 366 367 368 369 370 371
## 7.858032 9.853711 10.770694 11.860553 11.582459 13.304214 13.852817 14.336920
## 372 373 374 375 376 377 378 379
## 10.818551 10.468502 10.019509 9.677934 9.355059 9.031490 13.464253 11.100662
## 380 381 382 383 384 385 386 387
## 10.828725 10.901613 10.572265 10.469195 10.385394 9.574611 8.311332 7.455373
## 388 389 390 391 392 393 394 395
## 9.415499 8.811789 9.272832 19.516937 20.275463 19.452284 20.563896 23.256357
## 396 397 398 399 400 401 402 403
## 28.370741 30.397420 25.081901 27.482936 41.577907 29.373670 30.524598 32.153403
## 404 405
## 31.108659 34.366350Much lower values for archosaurs and dinosaurs – eggshell thickness seems to have become thinner for a given egg mass at the base of Ornithodira, and later increased in theropods. However, this pattern is dependent on a very small sample of pterosaurs, ornithischians, and sauropods; it is likely that the addition of new specimens attributed to any of these clades would considerably change that pattern.
We can see a strong increase in geckos and in eufalconimorphs – the latter having already been identified in Legendre and Clarke (2021).
References
- Legendre, L.J., Clarke, J.A., 2021. Shifts in eggshell thickness are related to changes in locomotor ecology in dinosaurs. Evolution 75, 1415–1430.
- Legendre, L.J., Rubilar-Rogers, D., Musser, G.M., Davis, S.N., Otero, R.A., Vargas, A.O., Clarke, J.A., 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature 583, 411–414.
- Norell, M.A., Wiemann, J., Fabbri, M., Yu, C., Marsicano, C.A., Moore-Nall, A., Varricchio, D.J., Pol, D., Zelenitsky, D.K., 2020. The first dinosaur egg was soft. Nature 583, 406–410.
- Stein, K., Prondvai, E., Huang, T., Baele, J.-M., Sander, P.M., Reisz, R., 2019. Structure and evolutionary implications of the earliest (Sinemurian, Early Jurassic) dinosaur eggs and eggshells. Scientific Reports 9, 4424.